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We introduce a new class of exchange-correlation potentials for a static and time-dependent Density 
Functional Theory of strongly correlated systems in 3D. The potentials are obtained via Dynamical 
Mean Field Theory and, for strong enough interactions, exhibit a discontinuity at half filling density, 
a signature of the Mott transition. For time-dependent perturbations, the dynamics is described in 
the adiabatic local density approximation. Results from the new scheme compare very favorably to 
exact ones in clusters. As an application, we study Bloch oscillations in the 3D Hubbard model. 
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Time-dependent quantum phenomena hold an impor- 
tant place in today's condensed matter research. A major 
theoretical challenge in this field is to describe strongly 
correlated systems out of equilibrium. 

In the last decade, Time-Dependent Density Func- 
tional Theory (TDDFT) has gained favor as a compu- 
tationally viable, in principle exact time-dependent de- 
scription of materials HJ [5]. The basic TDDFT vari- 
able is the one-particle density n and a key ingredient is 
the time-dependent exchange-correlation potential Vxc, 
embodying the complexities of the many-body problem. 
TDDFT applied to strongly correlated systems is in its 
beginnings. Describing these systems in equilibrium with 
static density functional theory (DFT) |3] is already a 
difhcult task gj. TDDFT retains these difficulties, but 
also adds another hurdle: Since time enters explicitly the 
formulation, Vxc depends on the history of n (memory ef- 
fects) [312]. 

In equilibrium, an effective ab-initio method to de- 
scribe strong correlations is the LDA-fDMFT [SJini, com- 
bining DFT in the local density approximation (LDA) 
with Dynamical Mean Field Theory (DMFT)[7j. DMFT, 
which treats correlations nonperturbatively via a local 
self-energy E [S], is also at the core of the DMFT+GW 
[51 110|. another ab-initio method, which deals with non- 
local correlations within the GW approximation [TT] . 
These DMFT-based methods rely on Green's function 
formulations, and the practical feasibility (in a foresee- 
able future) of a nonequilibrium generalization is not 
easy to assess, since Green's-function propagation scales 
quadratically [T^HTl] with the simulation time. 

TDDFT dynamics involves only one time variable. It 
would thus be useful to have exchange-correlation poten- 
tials suitable for strongly correlated systems. In equilib- 
rium, they could offer a better start for Green's function 
based ab-initio schemes. Out of equilibrium, they could 
be used for adiabatic LDA [THj dynamics via TDDFT 
and possibly be improved by including memory effects, 
absent in the adiabatic LDA. 

In this Letter we suggest a novel avenue to deal with 



strongly correlated systems in 3D and out of equilib- 
rium, by combining DMFT with TDDFT. For model 
strongly correlated systems in ID, exchange-correlation 
potentials for DFT were introduced [TB1[T7], and a Bethe- 
Ansatz-based LDA (BALD A) for v^c was proposed. Such 
Vxc^^'^^ was then used to introduce an adiabatic scheme 
for the TDDFT of the ID Hubbard model [H], and the 
spin dependent case was considered in [19] . However, 
some interesting effects due to correlations are specific to 
3D materials [3]. This requires new exchange-correlation 
potentials for [TD]DFT, that we propose here to obtain 
via Dynamical Mean Field Theory. 

We illustrate our method using the inhomogeneous 3D 
Hubbard model; as an initial application, we look at the 
Bloch oscillations in this model. Our main findings are: i) 
for the homogeneous 3D Hubbard model, above a critical 
interaction J7*^°", Vxc becomes a discontinuous function 
of n at half-filling (the size of the discontinuity increases 
at larger U:s): this is how the Mott metal-insulator tran- 
sition manifests itself in Vxc', n) the time-dependent den- 
sities from TDDFT-DMFT in the adiabatic LDA (here- 
after referred to as Aj^^f^pj,) compare very well with 
the exact ones in clusters; the agreement deteriorates 
for significantly non-adiabatic/strong perturbations; iii) 
^DA/FT gives a good description of the correlation in- 
duced beats in the Bloch oscillations, but no clear signa- 
tures of the damped regime, a fact most likely due to the 
lack of non-adiabaticity in our exchange-correlation po- 
tentials. While explicit for the 3D Hubbard model, our 
results also provide insight into the scope of TDDFT- 
DMFT for real strongly correlated materials. 

The Model- The time-dependent Hamiltonian for the 
Hubbard model is 

HiT)^-t^cl^Cj^ -h^U,ni-fh,i+^e,ni^-\-W{T), (1) 

where (ij) denotes nearest neighbor sites, a =tii 3.iid 
fiia = cl^Cia is the local density operator. We take t = 1 
as energy unit. The subscript i in the onsite energy and 



2 



repulsion term Ui allows for possible inhomogeneities. 
The external potential in time-dependent calculations is 
W{t) = X^icr Wi{T)nicr., with T bciug the time variable. 

Dynamical Mean Field Theory (DMFT).- In this study 
we neglect magnetic phases and allow only for paramag- 
netic solutions (n^ — n^), to focus on pure Mott physics. 
DMFT maps a Hubbard model on a simple cubic lattice 
onto a local problem representing one of the lattice sites 
(site 0) surrounded by a bath which describes the rest 
of the lattice [5]. In practice, one introduces auxiliary 
degrees of freedom to recover a Hamiltonian description 
of the local problem by identifying the site with the 
impurity site of an Anderson impurity model (AIM): 
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where Himp — U no-fnoi~ ^jluq , fi is the chemical potential, 
and the parameters Vi , e; are determined self-consistently. 
Self-consistency with the original lattice is obtained by 
requiring the impurity single-particle Green function (in 
Matsubara space POI) G{iuJn) to be identical to the local 
lattice Green function with identical self-energy E(iw„), 
i.e. 



G{iu)n) = I deD{e)[iuJn + M 
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where E(ia;„) is obtained via the local Dyson equation 
and D{e) is the non-interacting lattice density of states. 

We solved the self-consistent AIM using the exact- 
diagonalization (Lanczos) algorithm |21| . which trun- 
cates the number of auxiliary degrees of freedom to a 
finite, small number Ng. The results shown {Ng — 8) 
are converged against Ng. Once at self-consistency, the 
density n = X]cr('^o<T), the average double occupancy 
d — (nQ-fngi) (and thus the potential energy per lattice 
site V = Ud) are evaluated as averages on the impurity 
site (0) of the AIM. The total energy per site is given by 
Edmft = Kdmft + V, where Kdmft is given by [20] 
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DFT for the 3D Hubbard model - In a spin- independent 
DFT for the Hubbard model [22], the total energy is: 

Ey[n] = To W + EhIti] + E^dn] + ^ Vext{i)nt, (5) 



where v^j-t denotes the static external field (in the nota- 
tion of Eq. ([T|, Vext{i) = Ci)- In EQ- ^ = J2a'"'i<^^ 
while To[n] and Eh — \J2i^i''^i ^^'^^ respectively, the 
kinetic energy of the non-interacting system and the 
Hartree energy. 

We use an LDA for E^c and Vxc- v^cii) — v^dni), 
where E^c is obtained from the homogeneous 3D Hub- 
bard model, our reference system. We employed DMFT 



ing the ground state energy of the reference system. The 
DMFT impurity solver introduced some noise in the nu- 
merical solution, especially at low densities (n < 0.2) 
or close to half-filling (n < 1). Thus we first smoothed 
the data, and then performed a polynomial fitting for 
0.2 < n < 1. For n < 0.2, instead of DMFT, we used an 
analytic, asymptotically exact, form for the ground state 
energy of the 3D Hubbard model [23] . Including the sub- 
leading term at low n, one gets E^c = (Sttc^ — U) n'^/4:+ 
An^/'^ , where = u-^+-y ^^"^ scattering length for 
the model, 7 — 0.1263t~^ for a simple cubic lattice, and 
A is a fitting parameter. The piecewise analytical expres- 
sion for Exc was then differentiated to obtain z)*(n), the 
exchange-correlation potential for < n < 1. Due to 
electron- hole symmetry, in the entire density range [0, 2], 
we have = 0{1 - n)w,(n) - 0{n ~ l)i;,(2 - n), 

where 9 is the step function. 

Results for Exc and Vxc from DMFT are in Figs. [T^-b 
for several U values. On increasing U, a change of cur- 
vature occurs in Exc for n w 1 (Fig. [T^), and a cusp 
develops above a critical value [7^^°** sa 14. This induces 
a discontinuity in Vxc (Fig. ^p), a manifestation of the 
Mott-Hubbard metal-insulator transition in a DFT de- 
scription. Such behavior is quite different from that of 
Vxc in the ID Hubbard model, where the discontinuity 
occurs for any U > [T^]. In Figs, [it-d we present 
results for (present work) and ^/^^-^ [H] for 

U = 8 and 24. For U = 8, only wf/^^^ is discontinu- 
ous. The DMFT and BALDA exchange-correlation po- 
tentials are rather different from each other in the entire 



density range, with 



,,BALDA 



DMFT I 



reflecting the 



difference between ID and 3D reference systems. These 
features are generic for any U value. In particular, for 
U > ?7*^°**, the discontinuity in v^^^^^"^ is considerably 
larger than in w^*^^^, as seen in Fig. [iji. 

To ease the numerics, we slightly smoothed near n = 1 
the Vxc-s for U > U^'""- The thus obtained was 

used in our initial, ground state DFT-LDA calculations. 
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FIG. 1: (Color online) a-b): Exchange-correlation energies 
Exc and potentials v^c for the homogeneous 3D Hubbard 
model, for several values of the interaction U. c); DMFT 
(thick solid curve) vs ID BALDA results (thin solid curve) 
for Vxc when U = 8. d): same as c), but U = 24. 
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via the Kohn-Sham (KS) equations 



(6) 



where T = —tJ2{ij) a ^ia^j'y ^'^d ip^ is the K-th single par- 
ticle KS orbital, with rii = X^k*^"^ Iv'kC*)!^- The effective 
potential VKs{i) = Wa;c(«) + ^'ext(«), with vuii) = 

\Uini the Hartree potential, and Vxcii) = v^,^^'^{ni). 

TDDFT for the 3D Hubbard model. - To perform 
TDDFT real-time dynamics of the Hubbard model [18] , 
one propagates in time the KS orbitals ^Pk,{t) via the 
time-dependent Kohn-Sham equations: 



(7) 



to get the density ^^(t) X^^^ \VK{h''')\^. In general, 
VKsihT) = VH{i,T) + Vxc{hT) +Vext{i,T) dcpcnds non- 
locally on the density via Vxc- In the adiabatic LDA 
considered below, a local dependence in space and time 
is assumed: VxdhT) — >■ v^J^^ ^'^ {ni{T)) . 

Exact vs adiabatic-LDA dynamics. - In Fig. [2] we com- 
pare the A^^fprp densities to the exact ones, for a simple 
cubic cluster with 5"^ sites and open boundary conditions. 
We consider a highly inhomogeneous case, a single inter- 
acting impurity in the center i = 0: Ui = U5io in Eq.([T]), 
which should be a rather severe test for an adiabatic 
LDA based on a DMFT of the homogeneous 3D Hubbard 
model. We also set Ci — eo^io and Wi^r) — WQ{T)Sio. Due 
to symmetry, only 10 (Ngy) out of the total (per spin) 
125 one-particle eigenstates, those with nonzero ampli- 
tude at i = (active states), determine the static and 
time-dependent density at i = 0, making the size of the 
exact configuration space manageable [24] . 

We considered two time profiles for wo{t): gaus- 
sian and rectangle-shaped. In the latter, the ramping 
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FIG. 2: (Color online) Exact, TDDFT and Hartree density for 
U = 8 (a-c) and [/ = 24 (d-f) at the central site of a 5 x 5 x 5 
Anderson impurity cluster. In (c) and (f) the cluster has 
A^e = 70 electrons, and eo ~ —2.66 and eo — —4, respectively, 
to attain an initial (and the same) density close to half-filling. 
Otherwise Ne — 40 [23] and eo = 0. The time-dependent 
perturbation wo{t) acts always only at the impurity site, and 
its shape is shown below (d-f) (see main text for details). 



up/down of the pulse is faster. For [7 = 8, Fig. |2^-c, 
there is a very good agreement between Aj^^fpj, and the 
exact results when the perturbation is actually present. 
Afterwards, when Wq(t) has returned to zero, the agree- 
ment somewhat deteriorates for the faster perturbation, 
indicating the presence of non-adiabatic, non-local ef- 
fects in the exact dynamics. These however are not 
described by ^£)mft- Similar considerations apply for 
U = 24, Fig(2]i-f. Here, wf/^^^ exhibits a gap, and yet 
the ^£)MFT performs rather well at low filling (panels 
d-e). The performance worsens considerably closer to 
half- filling, panel f), for two reasons. 

The first concerns the ground state densities in Fig. 
[2] (i.e. at r = 0). When obtained with v^J^^'^, these 
agree very well with the exact ones, save for panel f). To 
understand why, we determined via reverse-engineering 
the exact ground state Kohn-Sham potential v^^g for all 



cases a-f). For panels a-e), 



for i 7^ 0, i.e. 



is essentially local. For case f), we found large nonzero 
values of at i ^ 0. Such nonlocality is missing in our 
^ISa/ft- second reason is that in Fig. [2|'the Aj^^fprp 
density crosses the discontinuity in the Vxc at half-filling. 
Such discontinuity was determined from the infinite, ho- 
mogeneous 3D Hubbard model, but is expected to be 
significantly modified in an inhomogeneous small cluster: 
such change is missed by our This causes the 

disagreement in Fig. |2j^ (e.g., oscillations occur around 
no = 1 in the ^Ijmft densities, which are different from 
those in the exact curve). 

In Fig. [2]we also display results obtained with v^^^^^ 
and Vxc — (Hartree dynamics). While conceptually un- 
justified in 3D, the adiabatic BALDA is an helpful tool 
to asses the validity of ^fj^FT- Thereby we see that 
for [/ = 8 the adiabatic BALDA induces density oscil- 
lations due to the discontinuity at half-filling which, in 

^BALDA „,DMFT 
BALDA 



(but not in v. 



) exists for any U > 0. Also, 



is considerably stronger than the DMFT coun- 
terpart, see e.g. Fig. [2f: The adiabatic BALDA density 
does not cross no = 1, due to the too large discontinuity 



of 



BALDA 



This shows that a 



^Xc ■ -^"i" UxC 

eludes correlations in 3D is needed. 



which properly in- 
Finally, the Hartree 
dynamics is much worse than the Aj^^fprp one. Other 
densities and perturbations gave consistent results |25j . 

Block Oscillations. - As an application of Aji,^fprp, we 
discuss briefly the Bloch oscillations in the 3D Hubbard 
model. The Bloch oscillations are the response of par- 
ticles in a lattice to an electric field F [26]. With no 
interactions, in a lattice a uniform F linearly increases 
the particle momentum until Bragg reflections occur and 
an oscillatory current sets in with frequency cj — F, i.e. 
the Bloch oscillations. These have been observed in su- 
perlattices and, quite recently, in ultracold-atom systems 
|27j . The effect of strong interactions on the Bloch os- 
cillations has been theoretically investigated mostly for 
ID and 2D bosons, and for ID and inflnite-dimensional 
fermions \^M I28| : Depending on the Hamiltonian param- 
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efficient and accurate [TDJDFT treatments of the equi- 
librium and nonequilibrium behavior of 3D strongly cor- 
related materials. 

We thank C.-O. Almbladh and U. von Barth for discus- 
sions. C.V. is supported by ETSF(INFRA-2007-211956). 



FIG. 3: Beats regime, with U = 2, = = 8, F = 
10, kx = 1, kyz = 10. a) Bloch oscillations of the density cen- 
troid nc(r) b) Power spectrum of ric- The temporal patterns 
persist for the rest of the (larger than in figure) simulation 
interval. 

eters {U/t and F/t in Hubbard models) damped behav- 
ior, beatings, or highly irregular patterns can take place. 
Here, in analogy with cold-atom experimental setups, 
we consider a 3D Hubbard model with an asymmetric 
parabolic confinement. We use a 3D cluster of 33 x 5 x 5 
sites (along x, y, z) with open boundary conditions. To 
have a ground state atomic "cloud" elongated in the x di- 
rection, we take a potential = ^ \_kxx1 -\- kyziyf + zf)^ 
with kx < kyz- 

At T = 0^ , we start the time evolution with a perturba- 
tion Wi{T) ~ 9{t) [— -|- Fxi] . Our time-dependent 
indicator for the Bloch oscillations is the one-particle cen- 
troid nc(T) = ~ xo)n'i{T), with the particle 

density at site i, Np = '^^ni, and Xq the center of the 
cluster in the x-direction. We examined several regimes 
where beats and damping of the oscillations are expected 
(the inherent parameter values were chosen by adapting 
to 3D the corresponding ones for such regimes in ID). 
For large driving fields, we found beats in the A'^fpj, 
results (Fig. [3^), with a frequency u) k F and extra 
peaks with a splitting Ao; « U (Fig. |3|d). No unequiv- 
ocal signatures of damped oscillations were observed for 
several setups in the expected parameter regime: It is 
highly plausible that this is due to the lack of memory 
in the ^fj^/^^y, suggesting the need for improved, non- 
adiabatic exchange-correlations potentials. 

Conclusions. - We introduced a method to deter- 
mine new exchange-correlation potentials for inhomoge- 
neous strongly correlated systems in 3D. We used DMFT 
as "engine" for the many-body calculations, but one 
could equally well have considered other approaches, such 
as, e.g., the Gutzwiller approximation or the Quantum 
Monte Carlo method. For 3D clusters with Hubbard- 
type interactions, A^fprp and exact results are in good 
agreement. While, in general, both nonlocal and non- 
adiabatic effects are needed in Vxc, we have here shown 
that an Aj^j^jp^p^ i.e. a Vxc which correctly incorporates 
correlations in 3D, is an effective starting point for a 
TDDFT description of strongly correlated systems in 3D 
and out of equilibrium. Moreover, by involving only one 
time variable, our approach offers a clear computational 
advantage in describing the time evolution of large inho- 
mogeneous 3D systems. Hopefully, these attractive fea- 
tures will stimulate further studies and pave the way to 
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